Quantum Monte Carlo study of the two-dimensional fermion Hubbard Model 
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We report large scale determinant Quantum Monte Carlo calculations of the effective bandwidth, 
momentum distribution, and magnetic correlations of the square lattice fermion Hubbard Hamilto- 
nian at half-filling. The sharp Fermi surface of the non-interacting limit is significantly broadened 
by the electronic correlations, but retains signatures of the approach to the edges of the first Bril- 
louin zone as the density increases. Finite-size scaling of simulations on large lattices allows us to 
extract the interaction dependence of the antiferromagnetic order parameter, exhibiting its evolu- 
tion from weak-coupling to the strong-coupling Heisenberg limit. Our lattices provide improved 
resolution of the Green's function in momentum space, allowing a more quantitative comparison 
with time-of-fiight optical lattice experiments. 
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I. INTRODUCTION 

Originally introduced to explain magnetism and metal- 
insulator transitions in solids with strong electronic 
correlations and narrow energy bands, the underly- 
ing physics of the fermion Hubbard Hamiltonian^^^ re- 
mains a topic of considerable discussion. In two di- 
mensions, when the lattice is doped away from half- 
filling, do the fermions condense into a superconducting 
state? If so, what is the symmetry of the pairing order 
parameter?-'^""^^ Do charge inhomogeneities (stripes and 
checkerboards) emerge, and what is their interplay with 
magnetic and superconducting orders?^^"^^ 

In contrast to this uncertainty concerning the proper- 
ties of the doped lattice, the qualitative behavior at half- 
filling (one fermion per site) is much more well under- 
stood. The interaction strength U causes both the devel- 
opment of long range antiferromagnetic order (LRAFO) 
and insulating behavior. Even so, there are still some re- 
maining open fundamental questions, for example, in the 
precise way in which the model evolves from the weak- 
coupling to strong-coupling limits, especially in two di- 
mensions. 

At weak-coupling, one pictures the insulating behav- 
ior to arise from a Fermi-surface instability which drives 
LRAFO and a gap in the quasiparticle density of states. 
On the other hand, at strong coupling the insulating be- 
havior is caused by Mott physics and the suppression 
of electron mobility to avoid double occupancy. These 
points of view are clearly linked, however, since for large 
U /t the Hubbard Hamiltonian has well defined local mo- 
ments and maps onto the antiferromagnetic Heisenberg 
model with exchange constant J = At^ /U 

Developing an analytic theory which bridges these 
viewpoints quantitatively is problematic. Hartree Fock 
(HF) theory provides one simple point of view, but pre- 
dicts LRAFO at finite temperatures in two dimensions, 
in violation of the Mermin- Wagner theorem. In fact, 
even in higher dimension when the Neel tempertaure Tn 



can be nonzero, HF theory predicts Tn oc U instead of 
the correct ly: J = At^ /U . Sophisticated approaches 
such as the self-consistent renormalized theory,^°'^^ the 
fluctuation-exchange approximation,^^ and two-particle 
self-consistent thcory^'^ obey the Mermin- Wagner the- 
orem and provide a good description of the Hubbard 
Hamiltonian at weak- coupling, but fail for large U/t. A 
recent approach^^ based on the mapping to the nonlinear 
sigma modeP^ has made some progress in connecting the 
two regimes. 

The need to pin down the behavior of the two- 
dimensional (2D) half filled Hubbard model more quan- 
tititively, in a way which links the weak-coupling and 
strong-coupling limits, is particularly germane at present 
with the achievement of cooling and quantum degener- 
acy in ultracold gases of fermionic atoms. ^^"'^^ Such sys- 
tems offer the prospect of acting as an "optical lattice 
emulator" (OLE) of the fermion Hubbard model, allow- 
ing a precise comparison of experimental and theoreti- 
cal phase diagrams which is difficult in the solid state, 
where the (single band) Hubbard Hamiltonian provides 
only a rather approximate depiction of the full complex- 
ity of the atomic orbitals. Obviously, the achievement of 
this goal is one which requires accurate computations. A 
particular issue in the field of OLE concerns whether the 
temperature dependence of the double occupancy rate 
changes sign during the course of the evolution from weak 
to strong coupling. '^^ 

It is the intent of this paper to present considerably 
improved results for the effective bandwidth, momentum 
distribution, and magnetic correlations of the square lat- 
tice fermion Hubbard Hamiltonian. We will employ the 
determinant quantum Monte Carlo (DQMC) method, 
which provides an approximation-free solution of the 
model, on lattices large enough to use finite-size scaling 
to, for example, reliably extract the antiferromagnetic 
order parameter function of interaction strength. 

There is a considerable existing body of QMC studies of 
the two-dimensional half-filled Hubbard model, both on 
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FIG. 1; (Color online) [(a) and (b)] The momentum distri- 
bution, Eq. (2), is shown for interaction strengths U ranging 
from U = 2t (one quarter the bandwidth) to U = W — 8t. A 
sharp Fermi surface is seen at weak-coupling as the momen- 
tum cuts across the Fermi surface at k = (7r/2,7r/2). Larger 
U broadens n(k) considerably. The occupation becomes sub- 
stantial outside the nominal Fermi surface. Panel (a) shows 
the full BZ, while panel (b) provides higher resolution for 
the portion of the cut perpendicular to the Fermi surface at 
(7r/2, 7r/2). (c) At (7 = 2t and /3t = 32, n(k) has only a weak 
lattice size dependence, apart from the better resolution as 
L increases, (d) For U = 2t on a 20 x 20 lattice, n(k) is 
converged to its low temperature value once T < t/8. (By 
contrast, the spin correlations reach their ground state values 
only at considerably lower T.) 



finite lattices and in infinite dimension. A partial list 
includes Refs. 13, 18, and 33-40. 



II. MODEL AND COMPUTATIONAL 
METHODS 



The fermion Hubbard Hamiltonian, 



rif- J) ("ii-^j -mE"-' 



(1) 



describes a set of itinerant electrons, represented by 
Cjg.(cj^), the annihilation (creation) operators at lattice 
site J and spin a. The corresponding number operator 

Hi, 



Cj^Cj^. The first term represents the hopping (ki- 



netic energy) of the electrons. We will choose the param- 
eter t = I to set our unit of energy. The non-interacting 
bandwidth W = 8t. U is the on-site repulsion of spin-up 
and spin-down electrons occupying the same lattice site, 
and fi is the chemical potential which controls the parti- 
cle density. We will mostly be interested in the properties 



of the model on TV = L x L square lattices at half-filling 
(the number of particles is equal to the number of lat- 
tice sites) which occurs at /i = with our particle-hole 
symmetric choice of the representation of the interaction 
term. 

We will also focus exclusively on the case of the square 
lattice. This particular geometry has several interest- 
ing features. The half-filled square lattice Fermi surface 
exhibits perfect nesting, and the density of states is (log- 
arithmically) divergent. As a consequence, the antifer- 
romagnetic and insulating transitions occur immediately 
for any nonzero value of the interaction strength U, in- 
stead of requiring a finite degree of correlation, as is more 
generically the case. 

Our DQMC algorithm is based on Rcf. 41 and has 
been refined by including "global moves" to improve 
crgodicity"*^ and "delayed updating" of the fermion 
Green's function, ^'^ which increases the efficiency of the 
linear algebra. Details concerning this new code are avail- 
able at Ref. 44. Some other approaches to fermion Hub- 
bard model simulations are contained in Refs. 39 and 
45-47. 



III. SINGLE-PARTICLE PROPERTIES 



We begin by showing single-particle properties. The 
momentum distribution n(k) — ^ (•^ko-'^ko-) 



tained directly in DQMC via Fourier transform of the 
equal-time Green's function Gji = (c- C;^) 



n(k) = 1 



2N ^ 



,*k-(j- 



(c. J ) 



(2) 



At [/ = and at half- filling, n(k) = 1(0) inside (outside) 
a square with vertices (7r,0), (0, tt), (— tt, 0), and (0, — tt) 
within the Brillouin zone (BZ). In Fig. 1(a), we show 
7i(k) around the complete BZ, while Fig. 1(b) focuses on 
the region near the Fermi-surface point (tt, 2,7r/2). In- 
teractions broaden the J7 = Fermi surface considerably. 
Figure 1(c) shows that data for different lattice sizes fall 
on the same curve. Smearing due to finite temperature 
effects is seen in Fig. 1(d) to be small below T = t/ 

m = 8). 

Recent optical lattice experiments^^ have imaged this 
Fermi surface for a three-dimensional cloud of fermionic 
■^'^K atoms prepared in a balanced mixture of two hy- 
perfinc states which act as the Hubbard Hamiltonian 
spin degree of freedom. In Fig. 2 we show a sequence 
of color contour plots for different densities at weak and 
intermediate couplings, U/t = 2,4. As in the experi- 
ments, and in agreement with Fig. 1, the Fermi surface 
may still be clearly discerned, and evolves from a cir- 
cular topology at low densities into the rotated square 
as the BZ boundaries are approached. Because of the 
sign problem"*®'"*^ which occurs in the doped system, the 
temperatures shown in the figure are rather higher than 
those used in Fig. 1 at half-filling. Another single-particle 
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FIG. 2: (Color online) Color contour plot depiction of the momentum distribution 7i(k) and its gradient Vn(k). (a) Left to 
right, n(k) at weak-coupling U — 2t and fillings p — 0.23, 0.41, 0.61, 0.79, and 1.0. (b) Vn(k) for the same parameters, [(c) 
and (d)] Intermediate coupling U = 4t and fillings p — 0.21, 0.41, 0.59, 0.79, and 1.0. The increased breadth of the Fermi 
surface with interaction strength is evident. The lattice size = 24 x 24 and inverse temperature fit = 8 except at [/ = 4t and 
fillings p = 0.59 and 0.79, where the sign problem restricts the simulation to inverse temperatures /3t = 6 and 4, respectively. 
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FIG. 3: (Color online) As the interaction energy U increases, 
the kinetic energy declines. Here we show t^a/t, the ratio 
between the expectation value of (cj^^^Cj^) at U with its 
value at f/ = 0, for a 10 x 10 lattice. Strong-coupling and 
perturbative graphs are also shown for fit = 12. 



quantity of interest is the effective hopping, 



(3) 



which measures the ratio of the kinetic energy at finite U 
to its non-interacting value. As the electron correlations 
grow larger, hopping is increasingly inhibited, and tcs is 
diminished. In Fig. 3, we show a plot of this ratio as 
a function of U for a 10 x 10 lattice. Note that despite 
the insulating nature of the system, the effective hopping 
is nonzero and does not serve as an order parameter for 
the metal-insulator transition. Indeed, toff is responsible 
for the superexchange interaction which drives antiferro- 
magnetic order. The effective hopping can be evaluated 
analytically at small and large U^^. The DQMC data 
interpolates between these two limits. 



IV. MAGNETIC CORRELATIONS 



We turn now to two-particle properties, focusing on the 
magnetic behavior. The real-space spin-spin correlation 
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FIG. 4: (Color online) The local moment {m?) is the zero 
spatial separation value of the spin-spin correlation function 
C(0, 0). In the non-interacting limit {m?) — |. As the in- 
teraction energy U increases, (m^) approaches 1, indicating 
the complete absence of double occupancy and a well-formed 
moment on each site. 

function is defined as 

C{\) - ((rij+it - "j+u)(«jt - "ji)) (4) 

and measures the extent to which the z component of 
spin on site j aligns with that on a site a distance 1 away. 
Although defined in Eq. (4) using the z direction, C(l) 
is rotationally invariant and in fact, we measure all three 
components to monitor ergodicity in our simulations and 
average over all directions to provide an improved esti- 
mator for the magnetic properties. 

The local moment {ni^) = C(0,0) = {{nj-\ - n-^iY) 
is the zero separation value of the spin-spin correlation 



FIG. 6: (Color online) Comparison of the absolute value of 
the equal-time spin-spin correlation function (7(1) at U = 2t 
and pt = 32 for L X L lattices with L = 8, 12, 16, 20, and 24. 
The inset is a close up view of the long-range correlations. 



function. The singly occupied states 1 1 ) and | \. ) have 
{m?) = 1 while the empty and doubly occupied ones | 0) 
and I ti)i have (m^) = 0. In the non-interacting limit, at 
half-filling, each of the four possible site configurations is 
equally likely. Hence the average moment {m?') = i. 
The on-site repulsion U suppresses the doubly occupied 
configuration and hence also the empty one, if the total 
occupation is fixed at one fermion per site. Ultimately, 
charge fluctuations are completely eliminated, {m?) 1 
and the Hubbard model maps onto the spin-^ Heisenberg 
Hamiltonian. This is illustrated in Fig. 4 for a 10 x 10 
lattice. By the time U = W = ^t, the local moment has 
attained 90% of its full value. Thermal fluctuations also 
inhibit local moment formation but the data shown for 
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FIG. 5: (Color online) Comparison of the equal-time spin-spin 
correlation function (7(1) on a 20 x 20 lattice with (n) = 1 and 
U = 2t for inverse temperatures pt = 12, 20, and 32. The hor- 
izontal axis follows the triangular path on the lattice shown 
in the inset. Anti-ferromagnetic correlations are present for 
all temperatures, and saturation is visible at pt = 32. 



FIG. 7: (Color online) The equal-time spin-spin correlation 
function (7(1) on a 24x24 lattice with {n} — 1. Data are shown 
for various U at low temperatures. Antiferromagnetic corre- 
lations are enhanced for larger values of U. The increase in 
statistical fluctuations with interaction strength in the DQMC 
method is evident. 
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FIG. 8: (Color online) The magnetic correlation function ^(k) 
for U = 2t and pt = 32. The horizontal axis traces out the 
triangular path shown in the inset. The function is sharply 
peaked at k = Q = (vr, n). 



different temperatures in Fig. 4 indicate they are mostly 
eliminated by the time T decreases below t/12 — W/96. 

Local moments provide an intuitive picture of the onset 
of long-range correlation in the strong-coupling regime. 
They first form on the temperature scale U, which acts 
to eliminate double occupancy, and then, at yet lower 
T, they order via antiferromagnetic exchange interaction 
with J = At'^/U. In contrast with this situation, weak- 
coupling correlations are better described as arising from 
the instability of the Fermi g gainst formation of a 
spin-density wave, a peculiarity of the square lattice, sug- 
gesting an ordering temperature proportional to U. 

Figure 5 shows the spin-spin correlation function in the 
latter regime {U = 2) for a 20 x 20 lattice at I3t = 12, 20 
and 32. The correlations extend over the entire lattice 
even at /3t — 12, i.e., the correlation length has become 
comparable to the system size already at this tempera- 
ture. The values of C(l) continue to grow as T is in- 
creased further, saturating at pt ~ 32. This observation 
disproves the commonly held idea that on finite clusters, 
the order parameter stop growing after the correlation 
length exceeds the linear size of the system. Such satu- 
ration happens at a much lower temperature, only after 
thermal fluctuations have been largely eliminated. 

A comparison of |C(1)| for [/ = 2 and different lattice 
sizes is given in Fig. 6, where data for L = 8 up to L = 24 
are plotted and we have taken the absolute value to make 
the convergence with L clearer. We have fixed (3t = 32 so 
the spin correlations have reached their asymptotic low- 
temperature values. As expected, the smallest lattice 
sizes (8 X 8) overestimate the tendency to order, with 
C(l)| significantly larger than values for larger L. How- 
ever, by the time L = 20 the finite-size effects are small. 

We next compare the spin-spin correlation function 
for various U at low temperatures on a 24 x 24 lattice 



FIG. 9: (Color online) The antiferromagnetic structure factor 
5(Q) as a function of inverse temperature at [/ = 2t for Lx L 
lattices with L = 4, 6, 8, 10, 12, 14, 16, 18, and 20. 

in Fig. 7. Long-range order is present at all interaction 
strengths. For each U, wc have chosen temperatures such 
that the ground state has been reached for this lattice 
size. Since statistical fluctuations increase significantly 
with U and with /3 in DQMC, it is advantageous not to 
simulate unnecessarily cold systems. As discussed above, 
such temperature should increase with U in the weak- 
coupling regime and scale proportionally to 1/U in the 
strong-coupling one. We indeed find the highest satura- 
tion temperatures in the intermediate regime, at U/t ~ 4. 

The magnetic structure factor iS'(k) is the Fourier 
transform of the real-space spin-spin correlation function 

5(k) = 5^e^'^'C(l), (5) 
1 

were S'(k) is plotted in Fig. 8 as a function of k for several 
lattice sizes with U = 2t and pt = 32. S(k) is small and 
lattice size independent away from the ordering vector 
k = Q = (7r,7r). The sharp peak at Q emphasizes the 
antiferromagnetic nature of the correlations on a half- 
filled lattice. 

In order to understand the implications of the lattice 
size dependence at the ordering vector in Fig. 8, we show 
in Fig. 9 the antiferromagnetic structure factor for U = 2t 
as a function of inverse temperature for various L. As 
expected, as L increases, a larger value of /? is required 
to eliminate the low-lying spin-wave excitations and to 
saturate the structure factor to its ground state value. 

It is seen from Eq. (5) that S'(Q) will grow linearly with 
the number of sites N — if there is long-range antifer- 
romagnetic order. Huse^° has used spin-wave theory to 
work out the first correction to this scaling, 

L2 3 ^' 
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FIG. 10: (Color online) Scaling results for U = 5t and fit = 
20. (a) Convergence of the extrapolated value of S'n(Q) as a 
function of the size of the excluded cluster. See Eq. (9) for the 
definition oi S„. A large fraction of the bias is removed 

without loss in precision. C{L/2, L/2) is plagued by a much 
larger error bar. (b) Scaling of S'n(Q) for n = 0, 1, and 5 and 
C[L/2, L/2) as a function of the inverse linear lattice size. 
The extrapolation was performed via a linear least-squares fit 
in all cases. 



Here maf is the antifcrromagnetic order parameter, rriaf 
can also be extracted from the spin-spin correlation func- 
tion between the two most distant points on a lattice, 
C{L/2, L/2), with a similar spin-wave theory correction, 



C(L/2,X/2) = ^ 



(7) 



We expect that the correction b < a since the struc- 
ture factor includes spin correlations at short distances 
which markedly exceed m^f , in addition to the finite lat- 
tice effects at larger length scales. For similar reasons 
S is expected to show larger corrections to the asymp- 
totic 1/L scaling behavior than C. Part of the origin of 
these corrections is trivial and evident from Fig. 6: S 
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FIG. 11: (Color online) Order parameter mat as a function of 
the interaction strength U. Earlier DQMC values (Ref. 51) 
are circles and Hartree-Fock theory scaled by the Heisenberg 
result from Ref. 48 are shown as a line with long dashes. The 
line with short dashes and the dashed and dotted line indicate 
the strong-coupling Heisenberg limits from QMC (Ref. 52 and 
spin-only low-energy theory with ring exchange (Ref. 53), re- 
spectively. Also shown is the RPA calculation of Ref. 54 (solid 
line). 



is the average of quantities as different as C(0,0) and 
C(L/2, L/2). In the large L limit, S and its 1/L finite- 
size error are dominated by the contribution of the large 
distance correlations but for small L a bias roughly pro- 
portional to 



C(0,0) -C(L/2,L/2) 
L2 



(8) 



is clearly present. This bias is larger for small U since 
the numerator in Eq. (8) gets smaller with increasing 
U and saturates to the value of the Heisenberg model at 
U /t c^S. On the other hand, the error bars on S are often 
significantly smaller than on C, a fact that is certainly 
advantageous in the final finite-size scaling analysis. 

A measure of magnetic order that incorporates the ex- 
tended linearity of C and the better statistical property 
of S is given by 



5„(Q) = 



E 



(9) 



where n is the number of distances shorter than 1^. Equa- 
tion (9) is nothing but the interpolation between S, cor- 
responding to n = 0; and C , the case of n = L^ — 1. Fig- 
ure 10(a) shows how the L — > od extrapolation evolves 
by increasing l^. When Ic is small the linear extrapo- 
lation is significantly biased by the small-L results. As 
Ic increases one reaches statistical convergence already 
for Ic = 1 (corresponding to n — 5, the cluster formed 
by the origin and the nearest neighbors) with minimal 
loss of statistical precision. That = 1 is all is needed 
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to reach statistical convergence is also manifest in Fig. 6 
where C(l) drops to almost a constant beyond this value. 

In Fig. 10(b) we show S'„(Q)/i^ versus inverse linear 
lattice size 1/L for U = 5t and n = 0, 1, and 5 since as 
shown in Fig. 10(a), there is no real gain in accuracy by 
excluding larger subclusters. The inverse temperatures 
is (3t = 20, so that S has reached its zero-temperature 
value regardless of L. We have repeated this finite-size 
scaling analysis for couplings U/t = 2,3,4,6,7, and 8 and 
extrapolated to infinite L using a linear least-square fit 
in 1/L. In Fig. 11, we show the resulting antiferromag- 
netic order parameter TOaf as a function of U /t employing 
the same normalization convention of the other Hubbard 
model studies reported in this figure. The early DQMC 
values obtained by Hirsch and Tang,^^ which are consis- 
tent with the ones obtained here, are shown. 

Figure 11 also summarizes a number of the avail- 
able analytic treatments. The line with long dashes is 
the result of Hartree-Fock theory scaled by the Heisen- 
berg result at strong coupling.^® The solid line is the 
random-phase- approximation (RPA) treatment in which 
the single-particle propagators in the usual RPA sum 
are also dressed by the one-loop paramagnon correction 
to their self-energy.^^ Also shown (line with dots and 
dashes) are the results of a spin-only low energy theory^'^ 
which includes not only the usual Heisenberg J = At^ /U 
but also all higher order (e.g., ring exchange) terms up to 
f^/U^. Finally, the line with short dashes is the Heisen- 
berg value determined by Sandvik.^^ 

V. SUMMARY 

In this paper we have presented the results of determi- 
nant quantum Monte Carlo calculations for the magnetic 
properties of the half-filled square lattice Hubbard Hamil- 
tonian. DQMC allows us to bridge the weak-coupling and 
strong-coupling regimes with a single methodology and a 
particular outcome of our work has been the calculation 
of the antiferromagnctic order parameter in the ground 



state as a function of U /t. We expect these values will 
be useful in validating OLE experiments on the fermion 
Hubbard model. 

By using an improved DQMC code, we have been also 
able to provide results on larger lattices than those orig- 
inally explored. This not only has allowed us to do 
more accurate finite-size scaling for the order parameter 
but we also obtain considerably better momentum resolu- 
tion and hence a description of the Green's function G(k) 
which also offers the prospect of improved contact with 
time-of-flight images from optical lattice emulators. ^^"■^^ 

This study demonstrates a significantly improved capa- 
bility to simulate interacting fermion systems, driven by 
more powerful hardware as well as algorithmic advances. 
Systems of 500 sites (fermions) can now be handled on 
a modest cluster of desktop computers. Larger system 
simulations can easily be contemplated using more pow- 
erful hardware and would scale as the cube of the number 
of particles, in the absence of the sign problem. This re- 
maining sign problem bottleneck prevents the study of 
the densities of most interest to high-temperature super- 
conductivity, i.e., dopings of 5 — 15 % away from half- 
filling and motivates the interest in analog computation 
for the Hubbard Hamiltonian.^^ It should be noted, how- 
ever, that the sign problem can be rather modest for 
other densities, e.g., quarter filling, where we now have 
the capability to undertake large scale studies. 
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